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The two-dimensional surface of a bi-axial ellipsoid is characterized by the lengths of its major 
and minor axes. Longitude and latitude span an angular coordinate system across. We consider 
the egg-shaped surface of constant altitude above (or below) the ellipsoid surface, and compute the 
geodetic lines — lines of minimum Euclidean length — within this surface which connect two points 
of fixed coordinates. This addresses the common "inverse" problem of geodesies generalized to 
non-zero elevations. The system of differential equations which couples the two angular coordinates 
along the trajectory is reduced to a single integral, which is handled by Taylor expansion up to 
, fourth power in the eccentricity. 
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The common three parameters employed to relate Cartesian coordinates to an ellipsoidal surface are the angles of 
latitude and longitude in a grid on the surface, plus an altitude which is a shortest (perpendicular) distance to the 
surface. The well-known functional relations (coordinate transformations) are summarized in Section HT1 
i-pH ■ The inverse problem of geodesy is to find the line embedded in the ellipsoid surface which connects two fixed points 
subject to minimization of its length. We pose the equivalent problem for lines at constant altitude, as if one would 
ask for the shortest track of the center of a sphere of given radius which rolls on the ellipsoid surface and meets two 
points of the same, known altitude. In Section [III[ we formulate this in terms of the generic differential equations of 
geodesy parametrized by the Christoffel symbols. 

In the main Section IIV1 the coupled system of differential equations of latitude and longitude as a function of path 
length is reduced to one degree of freedom, here chosen to be the direction along the path at one of the fixed terminal 
points, measured in the topocentric horizontal system, and dubbed the launching angle. Closed-form expressions 
of this parameter in terms of the coordinates of the terminal points have not been found; instead, the results are 
■ presented as series expansions up to fourth power in the eccentricity. 

The standard treatment of this analysis is the projection on an auxiliary sphere; this technique is (almost) completely 
ignored but for the pragmatic aspect that the case of zero eccentricity is a suitable zeroth-order reference of series 
' expansions around small eccentricities. 

o \ 

' II. SPHEROIDAL COORDINATES 

X: 

H A. Surface 

The cross section of an ellipsis of equatorial radius p e > p p with eccentricity e in a Cartesian (x, z) system is: 

Pl = pl{l-e*); 4 + 4 = 1 - (1) 

Pe Pp 

The ellipsis defines a geocentric latitude </>' and a geodetic latitude (f>, the latter measured by intersection of the normal 
to the tangential plane with the equatorial plane (Figure [TJ. On the surface of the ellipsoid [HI, §IX]Q0, §140]: 



Z rp" z p^ 

tan<^= — |; - = tan <j>' = -| tan <f>. (2) 
x Pp x pj 
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FIG. 1: Major semi-axis p e , minor semi-axis p p , straight distance to the center of coordinates p, geocentric latitude 0', geodetic 
latitude 0, and their pointing difference v. 



Supposed p denotes the distance from the center of coordinates to the point on the geoid surface, the transformation 
from (x, z) to is 

,/ p e cos0 . , p e {l - e 2 )sin0 

x = pcos<t> = = , z = psin0 = =. (3) 

V 1 — sin <p y 1 — e z sm 

v = <f> — 0' defines the difference between the geodetic (astronomical) and the geocentric latitudes, 

e 2 sin(20) msin(20) . .„ .. . ,„ ,. .„ .. -> . ,„ ,. ,. ... 

t&nv = i- % — = —t—t = m sin 20 - m 2 sin (20) cos(20 + m 3 sin 20 cos 2 (20 + • • • ; (4 

2(1 -e 2 sin 0) l + mcos(20) 



where the expansion of the denominator has been given in terms of the geometric series of the parameter 



The Taylor series in powers of sin(20) and in powers of m are 



(5) 



m . m 2 (2 + m) . o m 2 (5 -I- 25m + 15m 2 + 3m 3 ) . 

■ sm(20) + — ^ -f sm 3 (20) + i — ^ '- sm J (20) + • • • (6) 



1 + m 6(1 + m) 3 40(1 + m) 5 

= sin(20)m - ^ sin(40)m 2 + - sin(60)m 3 - ~ sin(80)m 4 + - sin(lO0)m 5 H . (7) 

2 3 4 5 

A flattening factor / is also commonly defined [ill ]. 

Pe 

The reference values of the Earth ellipsoid adopted in the WGS84 [3] are 

/ = 1/298.257223563, p e = 6378137.0 m. (9) 

The numerical evaluation of ^ in terms of a Fourier series with this parametrization is, in units of radians and 
arcseconds: 

v = 0.0033584338 sin(20) - 0.56395388 x 10~ 5 sin(40) + 0.12626678 x 10~ 7 sin(60) + • ■ ■ (10) 
= 692."72669sin(20) - 1."16324 sin(40) + O."OO26Osin(60) + • • • . (11) 

Slightly different values would emerge according to IERS conventions of 2003 [HI] . 



FIG. 2: A point with Cartesian coordinates (x, y, z) has a distance z to the equatorial plane, a distance \J x 2 + y 2 to the polar 
axis, a distance h to the surface of the ellipsoid, and a distance N + h to the polar axis, measured along the local normal to 
the surface [12]. The vector points North at this point. 



B. General Altitude 



If one moves along the direction of <f> a distance h away from the surface of the geoid, the new coordinates relative 
to © are 0, i, E3, [3 US HI 

x = p cos <))' + h cos (j>; z — p sin <j>' + h sin <j), (12) 
which can be written in terms of a distance iV (</>), 

N(cb) = Pe (13) 

y 1 — e 2 sm q> 

as 

x = (N + h) cos <p; z= [N(l - e 2 ) + h}sm<j). (14) 

Rotation of (|14p around the polar axis with geographic longitude A defines the full 3D transformation between (x, y, z) 
and (A, <f>, h), 

[N(cj)) + h] cos (f> cos X \ 

[N(<p) + /i]cos<?!)sinA . (15) 
[N(<j))(l-e 2 ) + h] sin</> J 

The only theme of this paper is to generalize the geodetic lines of the literature 0, [E El; El; HU IH| to the case 
of finite altitude h =/= 0. The physics of gravimetric or potential theory is not involved, only the mathematics of 
the geometry. It should be noted that points with constant, non-zero h do not define a surface of an ellipsoid with 
effective semi-axes p e , p + h — otherwise the geodetic line could be deduced by mapping the problem onto an equivalent 
ellipsoidal surface Q. 



III. INVERSE PROBLEM OF GEODESY 

A. Topocentric Coordinate System 

A line of shortest distance at constant height h =const between two points 1 and 2 is defined by minimizing the 
Euclidean distance, the line integral 




4 



for some parametrization \(4>). The integrand is equivalent to a Lagrange function C(4>, X, dX/d<fi) or C{4>, A, d(f>/dX). 
The gradient of of (fT5")) with respect to A and </> defines the vectors e w that span the topocentric tangent plane 



- [N((f>) + ft] cos <f> sin A 
e\ = | W(<p) + ft] cos(/>cos A 




where a meridional radius of curvature [25 



(1 -e 2 sin^ 



- sin cos A [M(<^>) + ft] 
e, , = | — sin 4> sin A [M(0) + ft] 
cos [M (0) + ft] 



(17) 



)3/2 



1 — e 2 sin 2 . 



(18) 



is defined to simplify the notation. Building squares and dot products computes the three Gauss Fundamental 
parameters of the surface [13] 



e\ = E = (N + ft) 2 cos 2 , 
e A • e = F = 0; 

1-e 2 



G 



l 2 



N- 



1 — e 2 sin 



(M + ft) 2 



(19) 
(20) 

(21) 



Specializing to ft = wc get the You formulae 24] . E and G provide the principal curvatures along the meridian and 
azimuth [5j , and the coefficients of the metric tensor in the quadratic form of ds 2 , 



S = J ^EdX 2 + 2Fd\d<j) + Gd(j) 2 . 



(22) 



B. Christoffel Symbols 



Christoffel symbols are the connection coefficients between differentials de e of the topocentric axis and of the 
positions did 3 , in a generic definition 



de e = ^e a T%dx. . 



a,8 



This format is matched by first computing the derivative of e A with respect to A and 4> (at ft = const) : 



de x = 



and of es with respect to A and 



(M + ft) sin <\> sin A 
de^ = I — (M + ft) sin <j> cos A 




-(N + ft) cos 4> cos A 
- (AT + ft) cos sin A | dX + 




(M + ft) sin (j> sin A 
- (M + ft) sin 4> cos A 




dX 



V 



Ne 2 
Ne 2 



^--3Ne 2 (l-e 2 ) T 

CI (f> V ' ( 

*tt -3iVe 2 (l-e 2 ) 7 



1 — e 2 sin 2 0) 2 

sin 2 * \ — e 3 ) + /»] 



- (iV + ft) 



iVe 2 (l 



COS 



' (1 — e 2 sin 2 (p) 2 1 — e 2 sin 2 <f> 

The next step splits these two equations to the expanded version of (|23p , 

de x = (e A r^ + e^rf A )dA+(e A r^ + e^ A )# 

A<£ 



cos cos A \ 
cos 4> sin A 
sin d> 



de = (e x r^ + e rt )dA+(e A r^ + e r^)#. 



The eight T are extracted by evaluating dot products of the four vector coefficients in (l2~4"|) - ([2l)j) by e A and e^, 

0: 



1 AA 



1 <f>\ 



1 \4> 



4>4> 



(23) 



(24) 



(25) 

(26) 
(27) 

(28) 
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r AA = (N + h) cos (j) smcj) 



1 - e 2 sin 2 . 



r^ = r^ = -sin ( 



r«\ = 3iVe 2 (l - e 2 )sin</>cos< 



/i(l-e 2 sin 2 (f>) + N(l-e 2 ) 7 

h(l-e 2 sin 2 0) + 7V(l-e 2 )^ 
{N + h)(l-e 2 sin 2 0) cos^ ' 

1 

'[/i(l - e 2 sin 2 0) + N(l - e 2 )](l - e 2 sin 2 . 



(29) 
(30) 
(31) 



The Euler-Lagrange Differential Equations S J x \/ EdX 2 + GdtjP- = for a stationary Lagrange density C (at F = 0) 
become the differential equations of the geodesic [17j , in the generic format 



dV 
ds 2 



v dx^dx^ = 



/I// 



The explicit write-up 



d 2 A 

" 


x dX dX 
^ Txx ds~ds" 


A dXd(/) 
^ x *ds~ds~~ 


A dcf>dX 
^* x ds~Ts~ 




d^d^ _ Q ^ 
' ds ds 


d 2 4> 


j, dXdX 
b xx dsds~ 


dXd<p 
A * ds ds " 


d^dA 
<^ A ds ds " 




<fyd<p _ 
' ds ds 



simplifies with (|28jl to 



d 2 A , dA d</> _ 
d? + 21 ^ds"ds" = ° ; 



d4> 
ds 



0. 



(32) 

(33) 
(34) 

(35) 
(36) 



IV. REDUCTION OF THE DIFFERENTIAL EQUATIONS 
A. Separation of Angular Variables 

Decoupling of the two differential equations (|3"5"|) - (f3"6"]) starts with the separation of variables in 

_ _ 2r A # 

(IS 

Change of the integration variable on the right hand side from s to allows to use the underivative of (|3"0|) 

r A?i(^)# = ^slI^O) + h] cos</>} + const 



to generate a first integral 



Exponentiation yields 



where 



log ^ = — 2 log[(iV + /i) cos <ft\ + const, 
ds 



dX 



ds ° 3 (N + h) 2 cos 2 0' 



c 3 = ^- (JVi + /i) 2 cos 2 
ds |i 



(37) 



(38) 



(39) 



(40) 



(41) 



6 



is a constant for each geodesic. Is has been defined with the azimuth fa and iVi = N(fa) at the start of the line, 
but could as well be associated with any other point or the end point 2. C3 is positive for trajectories starting into 
eastwards direction, negative for the westwards heading, zero for routes to the poles. Moving the square of (|40f into 
([55]) yields 

d 2 <p sm(j) l-e 2 sin 2 2 3Ne 2 (l — e 2 ) sine/) cos 4> ( dfa 2 

(N + h) 3 cos 3 <j> h(l - e 2 sin 2 fa + N(l - e 2 ) ° 3+ [h{l - e 2 sin 2 fa + N(l - e 2 )](l - e 2 sin 2 fa \~cb , 



To solve this differential equation, we substitute the variable <p by its projection r onto the polar axis, 

t = sin0, (43) 

which implies transformations in the derivatives: 

dr ,d<f> . 

Ts =cos ^ (44) 



"r . , dd> dd> ,d T 

d? = - Sm ^dJ + C ° S ^ ; (45) 



cos0 d^ -d^ + T {d^) ~d^ + c~o~^\te) ~d^ + T^ 2 ~{ds~) ■ (46) 

We multiply (|4"2")1 by cos fa then replace d<p/ds and d 2 fa l ds 2 as noted above, 

drV, r l-e 2 r 2 37Ve 2 (l-e 2 )r f dr\ 2 _ 



ds 2 1-r 2 (iV + M 3 (l-7- 2 ) /i(l-e 2 r 2 )+iV(l-e 2 ) 3 - e 2 r 2 ) + JV(1 - e 2 )](l - e 2 r 2 ) ^ds 

d 2 r t - e 2 r 2 ) 2 + JV(1 - 4e 2 r 2 + 2e 2 + 4e 4 r 2 - 3e 4 ) / dr 



ds 2 1-r 2 [/i(l-e 2 r 2 )+iV(l-e 2 )][l-e 2 T 2 ] V ds , 

r 1 - e 2 r 2 



(JV + h) 3 (l - r 2 ) - e 2 r 2 ) + N(l - e 2 ) 3 
This is a differential equation with no explicit appearance of the independent variable s, 

(l_ T 2 )^!l + h^-e 2 r 2 ) 2 + N(l-e 2 )(l + 3e 2 ~'ie 2 T 2 )_fdr\ 2 , rc 2 1 - e 2 r 2 



c 2 = 0. (47) 



ds 2 [/i(l - e 2 r 2 ) + JV(1 - e 2 )] [1 - e 2 r 2 ] V ds J (N + /i) 3 h(l ~ e 2 r 2 ) + N(l - e 2 ) 

and the standard way of progressing is the substitution 



dr d 2 r dp 

ds ' ds 2 ^ dr 



Pi —=p^z> ( 48 ) 



_ dp fe(l- e 2 T 2 ) 2 + iV(l-e 2 )(l + 3e 2 -4e 2 T 2 ) rc 2 1 - ^ = „ f49) 

1 ^dr [fc(l-e 2 r 2 ) + 7V(l-e 2 )][l-e 2 r 2 ] ^ (N + hf h(l - e 2 r 2 ) + N{1 - e 2 ) ' V 1 

This is transformed to a linear differential equation by the further substitution P = p 2 , dP/dr = 2pdp/dr, 

l n _ 2 dP fe(l-e 2 r 2 ) 2 + ^(l-e 2 )(l + 3e 2 -4e 2 r 2 ) rc 2 1 - e 2 r 2 

2 l TJ dr [/i(l-e 2 r 2 )+iV(l-e 2 )][l-e 2 r 2 ] T (iV + h) 3 - e 2 r 2 ) + N(l - e 2 ) { ' 

The standard approach is to solve the homogeneous differential equation first, 

dP _ h(l - e 2 T 2 ) 2 + N(l - e 2 )(l + 3e 2 - 4e 2 r 2 ) 

dr~ T [l-T 2 ][/i(l-e 2 r 2 ) + 7V(l~e 2 )][f-e 2 r 2 ] ' ( ' 
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After division through P, the left hand side is easily integrated, and the right hand side (incompletely) decomposed 
into partial fractions, 

, „ „ f h(l-e 2 T 2 ) 2 + N(l-e 2 ){l + 3e 2 -4e 2 T 2 ) J 

lOgP = —2 1 T— „, r . ,. TTTZ ^TTTZ T^^dT 



[1 - T 2 }[h(l - e 2 T 2 ) + JV(1 - e 2 )][l - e 2 



-2r 



1-t 



rdr + 3 



-2e 2 r 
1 - e 2 r 2 



dr + 6he< 



h(l-e 2 T 2 )+N(l-e 2 ) 



dr 



log(l - t 2 ) + 3 log(l - e 2 r 2 ) - 2 log h(l - e 2 T 2 fl 2 + p e (l - e 2 ) 



const. 



(52) 



P = const ■ 



(1 -r 2 )(l-e 2 r 2 ) 3 



= const ■ 



(l-r 2 )(l-e 2 T 2 ) 2 



const ■ (1 — t 2 ) 



[^(l-e 2 r 2 )3/ 2 + p e (l-e 2 )] z [/i(l-e 2 r 2 ) + iV(l-e 2 )] 2 G(r) 

Solution of the inhomogeneous differential equation (|50[) proceeds with the variation of the constant, the ansatz 

(l-r 2 )(l-e 2 r 2 ) 2 



P = c(r) 



]2 ' 



[/i(l - e 2 r 2 ) + JV(1 - e 2 )] z 
Back insertion into (|5"0|) leads to a first order differential equation for c(r) , 

dc(r) 1 - r 2 2rc§ f - e 2 r 2 



(53) 



dr G(r) (f - t 2 )(N + h) 3 h(l - c 2 t 2 ) + N(l - e 2 ) ! 
which is decomposed into partial fractions 



dc(r) 
dr 



eN 



2er 



1 



-It 



1 



1 - e 2 r 2 (JV + hf 1 - r 2 (JV + hf 



1 



1-T 2 



The ensuing integral over dr is solved by aid of the substitution t 2 = u, 

1 



c(r) 



const. 



'{N + h) 2 {l-T 2 ) 

Back into (|53p — using const to indicate placement of any member of an anonymous bag of constants of integration, 



P = 



const — 



1 



C5 



(N + h) 2 {l -t 2 ) 
(l-r 2 )(l-e 2 r 2 ) 2 



-;(l-r 2 )(l-e 2 r 2 ) 5 



- e 2 T 2 ) + N(l - e 

,2 
3 



c 2 (l-e 2 r 2 ) 2 



2 ) + N(l - e 2 )] z (N + h) 2 [h(l - e 2 T 2 ) + N(l - e 2 )Y 



1 



c 5 - 



l-T 2 - 



(h + M) 2 \ (N+ h) 2 J " u G(r 

The subscript 1 denotes values at the starting point of the curve, 



1-T- I 
C5 — ; I 1 - 



E(t) 



Pi = CB 



cos 0i (1 — e 2 sin 0i) 



c 2 (f - e 2 sin 2 fii) 



C5 



- e 2 sin 2 0i) + JVi(l - e 2 )] (N t + h) 2 [h(l - e 2 sin 2 ^) + JVi(l - e 2 )] 
cos 2 0i(l - e 2 sin 2 0i) 2 (dA/ds) 2 (iVi + hf cos 4 0i(l - e 2 sin 2 (/H) 2 



2 =Pl 



[h(l -e 2 sin 2 



1 )+iV 1 (l-e 2 )]' 



[Mi 



9 ■ 2 

e z sm 



1 )+AT 1 (l-e 2 )]' 



(54) 
(55) 
(56) 

(57) 
(58) 



Solving for C5 yields 



C5 



dA 
ds |i 



cos 2 </>i(JVx + h) 2 + 



p 2 [h(l-e 2 sm 2 cj )1 ) + N 1 (l-e 2 )] 2 _ ( d\ 



cos 2 4>i ( 1 — e 2 sin 0i ) 2 



ds |] 



cos 2 0i(iVi + /i) 2 + 



Pi 



COS z (pi 



-Gi 



v ds|i/ 1 (iVi + /i) 2 cos 2 0i 



(59) 
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FIG. 3: Two examples of trajectories of dr/ds [left, equation (|61[) ] ° r ds/dr [right] as a function of r. They may or may 
not pass through a zero r m of (|61[) while connecting the starting abscissa n with the final abscissa T2- There are four basic 
topologies, depending on whether r m is positive or negative, and depending on whether the sign change in dr/ds is from + to 
— or from — to + at this point. 



Compared with the differential version of ([22]) . 

ds 2 = Ed\ 2 + Gd ( f> 2 ; 1 = e(^\ +G(^j , (60) 

we must have C5 = 1. Note that ([56]) is essentially a write-up for p 2 ~ {d<p/ds) 2 and could be derived quickly by 
inserting (|40| directly into (|60| . 
The square root of (|56p is 



= + LifT 1 _ 2 _ g = + v 1 ~ T = rfr 

P /i(l-e 2 T 2 )+^(l-e 2 )y T (7V + /i) 2 ft + M(r) ds' 1 ' 

The sign in front of the square root is to be chosen positive for pieces of the trajectory with dr/ds > (northern 
direction), negative where dr/ds < (southern direction). The value in the square root may run through zero within 
one curve, dr/ds = at one point r m , such that the square root switches sign there (Figure[3]). This happens whenever 

{l-r 2 m )[N{r m ) + h] 2 =cl (62) 

yields a vanishing discriminant of the square root for some (A, cj>), that is, whenever the difference A2 — Ai is sufficiently 
large to create a point of minimum polar distance along the trajectory that is not one of the terminal points. 



B. Launching Direction 

So far we have written the bundle of all geodesies through (Ai,0i) in the format (joTj) , which specifies the change 
of latitude as a function of distance traveled. The direction at point 1 in the topocentric system of coordinates is 
represented by C3. To address the inverse problem of geodesy, that is to pick the particular geodesies which also 
runs through the terminal point (A2, 4>2) m fulfillment of the Dirichlct boundary conditions, the associated change in 
longitude, some form of ((40]) . 

dA = 1 £3 

ds C3 (N + h) 2 (l-r 2 ) E' 1 ' 

must obviously get involved. The strategy is to write down r(A) or A(r) with C3 as a parameter, then to adjust C3 to 
ensure with what would be called a shooting method that starting at point 1 at that angle eventually passes through 
point 2. Coupling of A to r is done by division of ([BTj) and ([6"5]) , 



dr _dr ds _dr dX _ (N + h) (l-r )yjl-r 2 - (N _f hy2 

d\~ ds~d\~ ds' ds _± c 3 (h + M) ' *• ' 
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Separation of both variables yields 

± 



dr- 



c 3 (h + M) 



{N+h) 2 (l-T 2 )Jl 



(65) 



(N+h) 2 



For short paths, small |A 2 — Ai|, the integral is simply If one of the cases occurs, where the difference in A 
is too large for a solution, the additional path through the singularity r m is to be used [with some local cxtremum 
in the graph r(s)], and the integral is to be interpreted as f r m — f^ 2 . Taking the sign change and the symmetry 
of the integrand into account, this is J^ 2 +2 J^j" = J^" 1 + f^™, twice the underivative at r m minus the sum of the 
underivative at t\ and t-i . In the right plot of Figure [31 the difference is in including or not including the loudspeaker 
shaped area between ti and r m . 

r m is the solution of (162]) . the positive value of the solution taken xidr/ds > at 0i, the negative value if dr/ds < 
at 4>i. The squared zero of dr/ds, r^, is a root of the fourth-order polynomial which emerges by rewriting (|62D . 



e 4 h\r 2 m f + 2e 2 h 2 [p\ - h 2 + e 2 (c 2 - h 2 )] (r 2 J 



( p 2 _ tfy + 2e 2 ( 2h 4 _ 2c 2 h 2 _ 2f) 2 h 2 



2 \2 



^ 2 ' - 2 <~ 2 ' h 2 )+e 2 (-(c 2 -h 2 ) 2 +pl(4 + h 2 ) 



2 



[- (p 2 e -hY+4 (p 2 e + h 2 )+e 2 (- (c 2 -h 2 

+ [(p e + h) 2 - c 2 ] [(p e - h) 2 - 4] = 0. (66) 
Alternatively, r m admits a Taylor expansion by expanding the zero of (|62p in orders of e: 

2 



± T m = \ 1 — 



H- 



4Pe 

2H 3 



1 



34Pe 



8H S 



111 
H 



H 1 



- 1 



l-^ + 0(e% 



where some maximum distance 



H = p e 



(67) 



(68) 



to the polar axis has been defined to condense the notation. 
The integrand in (|65[) has a Taylor expansion in e, 



± / dr 



cz/{pe + h) 



C 3 p e [(Pe+h) 2 (2 



2c 2 ] 



(l-r 2 )Jl 



2( Pe + hy ( 



We integrate the left hand side of (|6"9")1 separately for each power of e up to 0(e 4 ), 



3/2 



0(e 4 



Am. 



(69) 



C3Pe 



arctan— % 

VT 2iJ 2 

czPr / r 



T 



arctan — — | e" 



16H 7 VT 3 /2 



(j^/2 [~ 2H2 4Pe + 4H 2 4peT + 24p e - 6i/ 3 c^T + 6H 5 T - 9H 5 T 2 - 4H 4 p e T + 6H 4 p e T 2 



H 2 [3H 3 - 2p e H 2 



3c 2 iJ + 4c 2 p e ] arctan —== ] c 



0(e 6 



where 



T = 1 - t z 



±(A 2 -Ai), (70) 



(71) 



is some convenient definition of the trajectory's distance to its solstice r m . The first term is not the principal value 
of the arc tangent but its steadily defined extension through the entire interval of r- values, \t\ < r rn . It is an odd 
function of r, and phase jumps are corrected as follows: the inclination at r = has (by inspection of the derivative 
above at r = 0) the sign of C3. So whenever the triple product sgn r sgn C3 sgn(arctan .) is negative, one must shift the 
branch of the arctan by adding multiples of 7rsgnr sgnc3 to the principal value. 

Whether this is to be taken between the limits T2 and t\ or as the sum of two components (see above) can be 
tested by integrating up to r m , J r m , where the value of the underivative at r m is given by (1/2) arctan 0, effectively 
j sgnr m sgnc3 after selecting the branch of the inverse trigonometric function as described above. 
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Equation (|70f is solved numerically, where C3/H is the unknown, where T\ t 2 and Ai^ are known from the coor- 
dinates of the two points that define the boundary value problem, and where e and H are constant parameters. 
The complexity of the equation suggests use of a Newton algorithm to search for the zero, starting from (|A19|) . 
C3/H ps cos</>i cos (^2 sin(A2 — \\)/ sinZ, as the initial estimate. 

An alternative is to insert the power series 



C3 



.(0) 



4V 



4V 



(72) 



right into (|65p . integrate the orders of e term- by-term, and to obtain the coefficients 

with the equivalent powers of the right hand side. c[ 

(2) 

recursively calculated from linear equations. c 3 , for example, is determined via 



■ , , Cg 4 -* etc. by comparison 
^ is given by (| A19|l : the coefficients of the higher powers are 



(0) 
H 2 



\ 



1-r 2 - 




1 - 



„(or 
^ 

H 2 



arctan ■ 



Ht 



1-t 2 - 



H 



(0) 

H 2 



joy 



H 2 



c (2) 
H 



= 0. (73) 



The corresponding equation for is already too lengthy to be reproduced here, so no real advantage remains in 
comparison with solving the non-linear equation (|69p . 

Once the parameter C3 is known for a particular set of terminal coordinates, \(<fi) is given by replacing T2 and A2 in 
(|70p by any other generic pair of values. Two other variables of interest along the curve, the direction and integrated 
distance from the starting point, are then accessible with methods summarized in the next, final two sub-chapters. 



C. Nautical course 



The nautical course at any point of the trajectory </>(A) is the angle k in the topocentric coordinate system spanned 
by (direction North) and e A (direction East), measured North over East, 

dr e e A e A e 

— = cosk— — + sm/i- — r = sinK— 1= + cosk— 74 
ds |eJ e A y/E VG 



dr/ds is the differential of (fT"5|) with respect to s, where d/ds — (dX/ds)(d/d\) + (d(j>/ds)(d/d(p), 

dr dX dd> 

— oc — e A + — e . (75) 
ds ds ds 

The sign oc indicates that the left hand side of this equation is a vector normalized to unity, but not the right hand 
side. 

sinn/yE dX/ds dX 
cos k/VG d(j)/ds dej)' 

Insertion of HI]), flU) and flUD yield 

M + h 1 1 dr/d(/> cos0 

d& dcf> dr dr dr_ 

dX dr dX dX dX 



(N + h) cos (j> K d * dT &l ' 



Mixing (I64| into this yields the course, supposed <f> and C3 are known 

C3 



tan«; = ± ° 3 (78) 
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D. Distance To Terminal Points 



An implicit write-up for the distance from the origin S\ measured along the curve is given by separating variables 
in EH): 



h(l-e 2 r 2 )+N(l-e 2 ) 



s = ± I dr- 

l T1 (l_ e 2 T 2)^ 1 _ T 2_ 



(79) 



(N+h) 2 



Power series expansion of the integrand in powers of e, then integration term-by-term, generate a Taylor series of 
the form 

v ' iTi 

In the notation (|71|) . the components of the underivative read: 



s*- -* = H arctan 



T T 

— , = H arctan — = ; 

v/l - r 2 - c 2 /H 2 VT 



(80) 



5 (2) = _Pe 

4 



(1 + cl/H 2 ) arctan -= + r 



3(l-r 2 )-c 2 /H 2 



T 



T 



(81) 



s(4) = £{ - ^/ H ? - m<%/H)*(p e /H) + 4(c 3 /H) 2 (p e /H) 3] arctan-^ 

+ ^ [2AT(c 3 /H) 4 - 2AT(c 3 /H) 2 + 63T 2 (c 3 / H) 2 - 8(p e / H)(c 3 / H) 6 + 8(p e / H){c 3 / H) 4 

- S( Pe /H)T(c 3 /H) 4 + &(p e /H)T(c 3 /H) 2 - 12(p e /H)T 2 (c 3 /H) 2 - 27T 2 + 30T 3 



}. (82) 



V. SUMMARY 



Computation of the geodetic line within the iso-surface of constant altitude above the ellipsoid is of the same 
complexity as on its surface at zero altitude. Although the surface is no longer an ellipsoid, mathematics reaches an 
equivalent level of simplification at which one integral is commonly expanded in powers of the eccentricity or flattening 
factor. We have done this first for the parameter which provides a solution to the inverse problem, then for two of 
the basic functions, distance from the starting point and compass course. 



APPENDIX A: REFERENCE: SPHERICAL CASE 



The limit of vanishing eccentricity, e = 0, simplifies the curved trajectories to arcs of great circles, and presents an 
easily accessible first estimate of the series expansions in orders of e. The Cartesian coordinates of the two points to 
be connected then are 

(cos 4>i cos Ai \ / cos 4>2 cos A2 \ 

cos 0i sin Ai ; S2 = (p e + h)\ cos 02 sin A2 . (Al) 
sin 0i J y sin 02 / 

The angular separation Z is derived from the dot product si ■ 83 = |si 1 1«2 1 cosZ, 

cos Z — sin 0i sin 02 + cos 0i cos 02 cos(Ai — A2); < Z < it. (A2) 

Each point s on the great circle in between lies in the plane defined by the sphere center and the terminal points, and 
can therefore be written as a linear combination 

s = asi+/3s 2 , 0<a,(3. (A3) 
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The square must remain normalized to the squared radius 

s 2 = a 2 s\ + (3 2 s\ + 2a(3s 1 ■ s 2 = a 2 (p e + h) 2 + (3 2 (p e + hf + 2a(3(p e + h) 2 cos Z, (A4) 
which couples the two expansion coefficients via 

a 2 +/3 2 + 2a/3cosZ= 1. (A5) 
One reduction to a single parameter £ to enforce this condition is 

<■"< - p= ^ _ ■ 0<e<V2. (A6) 



y/l + sin(2£) cosZ yT + sin(2£) cos Z 

In summary, the point (|A3|) on the great circle of radius p e + h between si and s 2 has the Cartesian coordinates 



cos t / cos 0i cos Ai \ s ., /cos 02 cos A 2 

(p e + h) cos 0i sin Ai + (p e + h) = cos0 2 sinA 2 I (A7) 

s/1 + sin(2£) cos Z ^ J + sin(2£) cos Z ^ sin(?! , 2 

cos0(£) cos A(£) 

(Pe + fe) I cos 0(0 sin A(0 I . (A8) 



sin 0(0 



The z-component of this, 



. cos C sin 0i + sin ^ sin 2 

sm0(O = , , =— , (A9) 

yl + sm(2£) cosZ 

allows to convert £ into ip. No ambiguity with respect to the branch of the arcsin arises since — ir/2 < < ir/2. The 
ratio of the y and ^-components demonstrates the dependence of A on £, 

, ... cos £ cos 0i sin Ai + sin £ cos 2 sin A 2 , A mi 

tanA(£) = t t : — —^-t t r> ( Al ° 

cos 4 cos 0i cos Ai +sm( cos 2 cos A 2 

where the arctan branch is defined by considering separately the numerator and denominator under the restrictions 
that no signs are canceled. The azimuth angle a between the points at £ = and at £ follows from s-si = (p e +h) 2 cos a, 

cos £ / cos0icosAi\ gin £ / cos0 2 cosA 2 \\ / cos0icosAi 

cos a = I — : cos 0i sin Ai H ; cos 02 sin A2 ■ cos0isinAi 



s/1 + sin(2£) cos Z \ sin(j)i J y'l + sin ( 2 £) cos Z \ sin ^ J J \ sin0i 
sin £ cos Z + cos £ 



y'l + sin(2£) cosZ 



< a < Z. (All) 



The path length along the great circle perimeter is simply the radial distance p e + h to the center of coordinates times 
the azimuth angle a measured in radians, 

/-rrr , , . / sin £ cos Z + cos £ „ . , 

^2 _ ^ + fc y = ( pe + fa ) arccos ■> ^ < s < (p e + h)Z. (A12) 

y 1 + sin(2£) cos Z 

To calculate the direction in the e\- e^-plane at the starting point S\, we employ £ as the parameter that mediates 
between A and 0: 

dX 1 dX 1 dX df 1 dX , da . , . 
= = 2. — / (A13) 

ds p e + h da p e + h d£ da p e + h d£ g?£ 

To calculate dA/d£, use the derivative of (|A10[) with respect to £, 

1 dX — sin £ cos 0i sin A x + cos £ cos 02 SU1A2 

cos 2 A d£ cos £ cos 0i cos Ai + sin £ cos 02 cos A2 

, • \ \ ~sin£cos0i cosAi + cos £ cos 2 cos A 2 

- (cos £ cos 0i sin Ai + sin £ cos 2 sin A 2 ) — ■ — — r-^r . ( A14) 

(cos £ cos 0i cos Ai + sin £ cos 02 cos A2) 
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In particular at the starting point, where £ = 0, A = Ai and 



1 dX cos 02 sin A2 cos 2 cos A2 

2~r^e = 1 \ cos0ismAi 7 r—^. 

cos 2 Ai a£ ij cos^icosAi (cos 0i cos \\Y 

By multiplication with cos 0i cos 2 Ai 

cos 0i — = cos 02 sin A2 cos Ai — sin Ai cos 02 cos A2 = cos 02 sin(A2 — Ai). (A15) 
"s |i 

To calculate d<j/d£, we convert the cosine in (|All|) to the sine, 

r. ~ — sin £ sin Z , A1 „\ 

Sincr = v l— cosher = — , Alo) 
y/l + sin(2£) cos Z 

The derivative of this with respect to £ is 

da _ . „ 1 1 . ^ . „ 2cos(2£)cosZ 

cos a— — cos 4 sin z — : — - sm 4 sm Z- 



d£ v/l + sin(20 cosZ 2 (1 + sin(2£) cos Z) 3 / 2 ' 

in particular at the starting point, a = £ = 0, 

% = sinZ. (A17) 

«? |i 

Insert this derivative and (|A15j) back into (|A13|1 

d\ 1 cos02 sin(A 2 — Ai) 1 



(A18) 



ds 1 1 p e + h cos 0i sin Z ' 

to obtain the master parameter C3 for the spherical case with (|41[) . 



4 o) = (p e + fe) COS(/>lCOS02 t (A2 ' Al) ; (e = 0). (A19) 
sin ^ 



APPENDIX B: ZERO ALTITUDE 



In particular on the ellipsoid, at h = 0, (formal) solutions to some integrals exist. An underivative of (|79p in terms 
of Elliptic Integrals of the Second Kind is [|, 3.158.11] 0, 219.07], 



S = ±p e 



'I 



Pi- 



er 



E arcsin( -) \ C - e r 

sm Q 



1-e 2 c 2 
1 — e 2 r 2 p 2 



+ const: h = 0, 



where an auxiliary modular angle C, is defined from 

sin£ = e\ 



1-4/ ft 

1 - c 2 e 2 /p 2 ' 



The left hand side of (|65|) can be written as an Elliptic Integral of the Third Kind [1, 3.157.7], 



± 



C3 



1-e 2 



=n 



sin £ 



eT 



L ; arcsin(— — -) \ C + const — A; /i = 0. 
e 2 sm £ ' 



(Bl) 



(B2) 
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APPENDIX C: NOTATIONS 



C5 , C3 constants of integration 

E, Ex Gauss parameter eq. (fH))) and its value at curve origin (Aj., <pi, h) 

E(. \ .) Incomplete Elliptic Integral of the Second Kind in the Abramowitz-Stegun notation [1, §17] 

e eccentricity of the ellipsoid, eq. ([1]) 

e\ t d> horizontal topocentric coordinate vectors at (A, </), h) 

f flattening, eq. (0) 

F Gauss parameter, eq. (|2"0f 

G, Gi Gauss parameter, eq. (|2ip and its value at curve origin 

T\ Christoffel symbols in (A, <f>, h) coordinates 

h vertical distance to surface of ellipsoid 

H maximum distance to polar axis (in the equatorial plane), eq. (|68[) 

k nautical angle, North over East in the topocentric tangential plane 

A longitude 

M a radius of curvature on the ellipsoid surface, eq. p8[) 

N a radius of curvature on the ellipsoid surface, eq. (fT3j) 

p distance ellipsoid center to foot point on the surface 

p e ,p equatorial, polar radius of ellipsoid, eq. (Q]) 

s distance along geodetic line measured from curve origin 

s vector from ellipsoid center to point on geodetic line 

S length of curved geodetic trajectory; equals S2 at (A2, h) 

sgn sign function, ±1 or 

a azimuth along great circle, eq. (|A12[) 

T normalized distance from closest polar approach, eq (f7Tj) 

t, T\ sin (f> and its value at start of curve, eq. (|43|l 

V geodetic minus geocentric latitude 

<j) geodetic latitude 

4>' geocentric latitude 

II(.; . \ .) Incomplete Elliptic Integral of the Third Kind [1, §17] 

x,y,z Cartesian coordinates from ellipsoid center, eq. (fT5|) 

£ parametrization of great circle (spherical case), eq. (|A7|) 

Z cone angle of circular section (spherical case), eq. (IA2|) 



APPENDIX D: FEM IMPLEMENTATION 



The simplest numerical solution of the inverse problem without restriction on the eccentricity could be a finite- 
element (FEM) integration of (p5|) and iterative adjustment of the single free parameter C3. The following Java 
program implements this approach. Member functions in overview: the constructors Geod define a surface from the 
parameters p e , e and h in which the geodetic line is embedded. getCartesian computes the vector (j 15[) . curvN 
computes (|13p . dNdtau computes 



f latt computes ((HJ. curvM computes ([75]). dMdtau computes 

dM ,»/M 
— = 3Mt- 

ar 1 — e z T z 



' /U SMW^. (D2) 



d2Mdtau2 computes 



GaussE computes (|19p . dEdtau computes 



d 2 M , N e 2 (l + 4e 2 r 2 ) , , 



^- = -2r(N + h)(M + h). (D4) 
ctr 
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d2Edtau2 computes the next higher derivative, d 2 E/dr 2 . GaussG computes (|2Tj) . tauTropic solves (j66]) for r m with a 
first estimate taken from (|67[) : the function is not actually called in this implementation, dtaudlambda computes dr/dX 
via (|64|) . referencing one factor to (|6Tj) . dsdlambdais ds/dX, the inverse of (|63|) . discrT computes 1 — r 2 — c|/(A^ + /i) 2 , 
which generalizes (|7ip to nonzero e. Its derivatives 



d 

d7 
£_ 

dT 2 



[N(T) + hf 



[N(t) + hf 



2t 



Ne 2 c 2 3 



(N + hf(l-e 2 T 2 ) 
1 



Nc 2 e 2 



(N + /i) 3 (l-e 2 r 2 ) 



3fteV 



2^2 



(iV + /i)(l-e 2 r 2 ) 



(D5) 
(D6) 



are implemented in dTdtau and d2Tdtau2. dtauds calculates (|6T|) . d2taudlambda2 calculates the derivative of (l64|) . 



dV 
dA 2 



d E(T)y/l — T 2 — c 2 /(7V + h) 2 dr d E(r)y/T 



c 2 /{N + hf 



dX 



c 3 (h + M) 



dX dr 



c 3 (h + M) 



(D7) 



d3taudlambda3 is the next higher order application of Bruno di Faa's formula to relegate derivatives d/dX to derivatives 
d/dr, [1, 0.430], 



dX 3 dX 2 dr 



ds2dlambda2 calculates 



c 3 {h+M) 



d 2 ^ 
dX 2 



d 3 r _ d 2 r d E(t)^/1-t 2 -c 2 /(N + h) 2 (dr\ 2 d 2 E{r)y/l 



4/{N + hf 



dX dr 2 



c 3 (h + M) 



d E 
dX c 3 



1 dr d 
c 3 dX dr 



-E. 



(D8) 



(D9) 



c3Sphere returns the estimate (IA19|) . dtauds Signum returns the sign of dr/ds at ti, obtained by considering the 
sign of the derivative of (|A9[) with respect to £. adjLambdaEnd modifies <p2 modulo 2ir to select the smallest value of 
1 02 — 01 1 - naut Angle computes k from (|78p . 

tauShoot walks along a geodetic line on a discrete mesh of width AA by extrapolating 



t\+ax ~ T\ 



dr 
dX 



AA 



l d 2 r 

2dy 



(AA) 2 



6dA3 



(AA) 3 , 



(D10) 



initialized at Ai,0i, given c 3 . The equivalent formula is used to build up sa+aa> currently only implemented up to 
the order (AA) 2 , since d3sdlambda3 returns 0. c3shoot calls tauShoot four times to adjust c 3 such that the error 
by which <f>2 was missed — returned by tauShoot — is minimized. The first call assumes (|Al9j) . the second takes an 
arbitrary small offset, and the third and fourth estimates are from linear and quadratic interpolations in the earlier 
calls to zoom into a root of this error as a function of c 3 . The last of these runs tabulates the Cartesian coordinates 
(fT5")) . A, 0, k and s on a subgrid of the A-mesh. main collects some adjustable parameters plus the pairs (Ai, <pi) and 
(X2,4>2), and calls c3shoot to solve the inverse problem of geodetics. 



package org.nevec . r jm 



/** Solution of the inverse problem of geodesy for biaxial ellipsoid. 

* The inverse problem of geodesy is solved given the parameters of 

* a biaxial ellipsoid (equatorial radius and eccentricity) , a common 

* altitude of a surface above the ellipsoid, and pairs of geodetic 

* coordinates (latitude, longitude) for a starting and a final point 

* of a trajectory. 

* <p> 
* 

* The first order differential equation of the smooth geodetic latitude 

* as a function of longitude is solved by moving from one point to the 

* next on a finite grid of equidistant sampling points on the longitude, 

* using a predictor-only Newton approximation. 

* The parameter which defines the initial direction in the local 

* tangential plane is gathered by a shooting method with a starting 

* guess obtained from a spherical approximation and a fixed number of 

* iteration through the parameter space. 
* 

* Ssince 2008-04-12 
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* @see preprint <a href ="http: //arxiv. org/abs/0711 . 0642">Geodetic line at constant altitude above the Ellipsoid</a> , 

* \protect\vrule widthOpt\protect\href {http : //arXiv . org/abs/0711 . 0642}{arXiv : 071 1 . 0642} [math .MG] 

* (Sauthor Richard J. Mathar 
*/ 

class Geod 
{ 

/** The inverse flattening factor of the IERS 2003 convention. 
*/ 

static public final double IERS_TN32_FLAT= 298.25642 ; 

/** The inverse flattening factor of the IERS 1996 convention. 
*/ 

static public final double IERS_TN21_FLAT= 298.25642 ; 

/** The inverse flattening factor of the GRS80 model. 
*/ 

static public final double GRS80_FLAT= 298.257222101 ; 

/** The inverse flattening factor of the WGS84 model. 
*/ 

static final double WGS84_FLAT= 298.257223563 ; 

/** Equatorial radius of the IERS 2003 convention, meters. 
*/ 

static public final double IERS_TN32_RH0_E=6378136 . 6 ; 

/** Equatorial radius of the IERS 1996 convention, meters. 
*/ 

static public final double IERS_TN21_RH0_E=6378136 . 49 ; 

/** Equatorial radius of the GRS 1980 system, meters. 
*/ 

static public final double GRS80_RH0_E=6378137 . ; 

/** Equatorial radius of the WGS84, meters. 
*/ 

static public final double WGS84_RH0_E=6378137 . ; 

/** Three enumeration values to convert angular units on user input and output 
*/ 

interface AngleUnit { 

int RAD =0, DEG = 1, G0N = 2; 

} ; 

/** Equatorial radius in meters. 
*/ 

double rho_equat ; 
/** Eccentricity. 

* For the Earth of the order of 0.08. 
*/ 

double eccen ; 

/** Altitude of the inverse problem above the ellipsoid. 

* Value equals on the ellipsoid surface. May have both signs. 

* Measured in the same units as Geod: :rho_equat . 
*/ 

double altit ; 

/** Default Constructor. 

* The variables are set to the defaults of the surface of the WGS84 geoid. 
*/ 

public GeodO 
{ 

this (WGS84_RH0_E , Math . sqrt ((2.0-1. /WGS84_FLAT) /WGS84_FLAT) , . ) ; 

} 

/** Constructor. 

* Oparam rho equatorial radius [m] 

* Oparam e eccentricity 



* Oparam h altitude above the reference ellipsoid [m] 

* Otodo catch error conditions of e>=l, rho<0 or h< polar radius. 
*/ 

public Geod(double rho, double e, double h) 
{ 

rho_equat = rho ; 
eccen = e ; 
altit = h ; 

} 

/** Convert geodetic to Cartesian coordinates. 

* Oparam tau the sine of the geodetic latitude 

* Oparam lambd the geodetic longitude [rad] 

* (Sreturn the three components of the Cartesian coordinates 

* in the same units as rho_equat and altit . 

* Osee eq (15) of the preprint 
*/ 

public double [] getCartesian(f inal double tau, final double lambd) 
{ 

double [] cart = new double [3] ; 
final double N = curvN(tau) ; 

final double sinphi = Math. sqrt (1 . -tau*tau) ; 
cart[0] = (N+altit)*sinphi*Math. cos (lambd) ; 
cart[l] = (N+altit)*sinphi*Math.sin(lambd) ; 
cart[2] = (N*(l.+eccen)*(l.-eccen)+altit)*tau ; 
return cart ; 

} 

/** Curvature parameter N [m] 

* Oparam tau sine of the latitude 

* Oreturn N according to eq (13) of the preprint. 
*/ 

public double curvN(final double tau) 
{ 

return rho_equat/Math. sqrt ( (1 . +eccen*tau) * (1 . -eccen*tau) ) ; 

} 

/** Derivative d N/d tau [m] . 

* Oparam tau sine of the latitude 

* Oreturn the derivative of N(tau) 
*/ 

double dNdtau(double tau) 
{ 

return curvN(tau) *eccen*eccen*tau/ ( (1 . +eccen*tau) * (1 . -eccen*tau) ) ; 

} 



/** Convert eccentricity to flattening factor. 

* Oparam e the (first) eccentricity 

* Oreturn f = 1-sqrt (l-eccentricity~2) . 
*/ 

public static double flatt (final double e) 
{ 

final double e2 = e*e ; 
if ( e2 < 3.8e-3) 

/* a Taylor approximation up to 0(e~14), relative accuracy to E-16. */ 

return e2* (0 . 5+e2* (1 . /8 . +e2* (1 . /16 . *e2* (5 . /128 . +e2* (7 . /256 . +e2* (21 . /1024. +e2*33 . /2048 .)))))) 

else 

return 1 . -Math. sqrt (1 . -e2) ; 

} 

/** Flattening factor. 

* Oreturn f = 1-sqrt (l-eccentricity~2) . 
*/ 

public double flatt () 
{ 

return flatt (eccen) ; 

} 
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/** Equatorial radius of curvature M [m] 

* (Spar am tau sine of the latitude 

* Oreturn M according to eq (18) of the preprint. 
*/ 

public double curvM(final double tau) 
{ 

return curvN(tau) * (1 . -eccen*eccen) / ( (1 .+eccen*tau) * (1 . -eccen*tau) ) ; 

} 

/** Derivative d M /d tau [m] 

* (Sparam tau sine of the latitude 

* Oreturn the derivative of curvMO with respect to the parameter tau. 
*/ 

public double dMdtau(final double tau) 
{ 

return 3 . *curvM(tau) *eccen*eccen*tau/( (1 .+eccen*tau) * (1 . -eccen*tau) ) ; 

} 

/** second Derivative d~2 M /d tau~2 [m] 

* Sparam tau sine of the latitude 

* Oreturn the second derivative of curvMO with respect to the parameter tau. 
*/ 

public double d2Mdtau2(f inal double tau) 
{ 

final double etau2 = eccen*eccen*tau*tau ; 

return 3 . *curvM (tau) *Math . pow (eccen/ ( 1 . -etau2) , 2 . ) * ( 1 . +4 . *etau2) ; 

} 

/** Gauss parameter E [m~2] . 

* Oparam tau sine of the latitude 

* Oreturn E according to eq (19) of the preprint. 
*/ 

double GaussE(final double tau) 
{ 

return Math.pow(curvN(tau)+altit , 2 . ) * (1 . -tau*tau) ; 

} 

/** Derivative d E/ d tau [m~2] 

* Oparam tau sine of the latitude 

* Oreturn the derivative of GaussEO with respect to the parameter tau 
*/ 

public double dEdtau(final double tau) 
{ 

return -2 . *tau* (curvN(tau)+altit) * (curvM(tau)+altit) ; 

} 



/** Derivative d~2 E/ d tau~2 [m~2] 

* Oparam tau sine of the latitude 

* Oreturn the 2nd derivative of GaussEO with respect to the parameter tau 
*/ 

public double d2Edtau2(f inal double tau) 
{ 

/* Use the Leibniz rule applied to the first derivative, which was -2*tau* (N+h) (M+h) 
*/ 

final double Nh = curvN(tau)+altit ; 
final double Mh = curvM(tau)+altit ; 

return -2 . * (Nh*Mh+tau*dNdtau(tau) *Mh+tau*Nh*dMdtau(tau) ) ; 

} 

/** Gauss parameter G [m~2] . 

* Oparam tau sine of the latitude 

* Oreturn G according to eq (21) of the preprint. 

* Osince 2008-04-22 
*/ 

double GaussG(final double tau) 
{ 

return Math.pow(curvM(tau)+altit , 2 . ) ; 

} 
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/** Minimum polar distance. 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn the value of tau_m according to eq (66) of the preprint that could be reached if the trajectory were 

* closed over the surface. That is the maximum of I sin (phi) I, or the value 

* of |sin(phi)| where the derivative of sin(phi) with respect to the longitude 

* becomes zero. 

* (Ssince 2008-04-22 
*/ 

double tauTropic (final double c3) 
{ 

final double c3sqr = c3*c3 ; 
final double H = rho_equat+altit ; 
final double rhoebar = rho_equat/H ; 
final double esqr = eccen*eccen ; 

/* Initial estimate from Taylor expansion in orders of squared eccentricity 

* See eq (67) of the preprint 
*/ 

double tau2 = (1 . -c3sqr/(H*H) ) * (l+c3sqr*rhoebar*esqr/(H*H)* 

(1 . +0 . 25* (3 . * ( 1 . -c3sqr/ (H*H) ) +rhoebar* (7 . *c3sqr/ (H*H) -3 . ) ) *esqr ) 
) ; 

final double hsqr = altit*altit ; 

final double rhosqr = rho_equat*rho_equat ; 

/* coefficients in the 4th order polynomial of tau~2 which is to be solved */ 
double [] coeff = new double [5] ; 

coef f [0] = (Math.pow(rho_equat+altit ,2 . ) -c3sqr) * (Math.pow(rho_equat-altit , 2 . )-c3sqr) ; 
coeff[l] = 2 . * (-Math.pow(rhosqr-hsqr ,2 . ) 
+c3sqr* (rhosqr+hsqr) 

+esqr* (-Math.pow(c3sqr-hsqr , 2 . )+rhosqr* (c3sqr+hsqr) ) 
) ; 

coeff [2] = Math.pow(rhosqr-hsqr ,2 . ) 

+2. *esqr*(2. *hsqr*hsqr-2. *c3sqr*hsqr-2 . *rhosqr*hsqr-rhosqr*c3sqr) 

+esqr*esqr*Math.pow(c3sqr-hsqr ,2 . ) ; 
coef f [3] = 2.* Math.pow(altit*eccen,2.)*(rhosqr-hsqr+esqr*(c3sqr-hsqr)) ; 
coeff [4] = Math.pow(altit*eccen,4. ) ; 

/* Small number of loops for self-consistency, using some Horner scheme for 

* the function to be zeroed and its first derivative */ 
for(int loop=0; loop < 4 ;loop++) 

{ 

System. out .println("# taum2 "+tau2) ; 

tau2 -= ( coeff [0]+tau2*(coeff [1] +tau2* (coef f [2] +tau2* (coef f [3] +tau2*coef f [4]))) 
)/( 

coeff [1] +tau2* (2 . *coef f [2] +tau2* (3 . *coef f [3] +tau2*4 . *coef f [4] ) ) 
) ; 

} 

return Math. sqrt (tau2) ; 

} 



/** Derivative d tau/d lambda. 

* Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn absolute value of the derivative according to eq (64) of the preprint. 
*/ 

double dtaudlambda(double tau, double c3) 
{ 

/* treat the factors (N+h) "2* (l-tau~2) and the sqrt ( . . ) /(h+M) separately */ 
return GaussE(tau) *dtauds (tau, c3) /c3 ; 

} 

/** Derivative d s/d lambda. 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn absolute value of the derivative of the path length according to eq (63) of the preprint. 
*/ 

double dsdlambda (double tau, double c3) 
{ 



return GaussE(tau) /c3 ; 

} 



/** Discriminant under the square root of d tau / ds . 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn l-tau~2-(c3/ (N+h) ) "2 . 
*/ 

double discrT(final double tau, final double c3) 
{ 

/* N+altitude */ 

final double Nh = curvN(tau)+altit ; 
return 1 . -tau*tau-Math.pow(c3/Nh, 2 . ) ; 

} 

/** Derivative of T with respect to tau. 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* Oreturn derivative of discrTO with respect to tau 
*/ 

double dTdtau(final double tau, final double c3) 
{ 

final double N = curvN(tau) ; 
final double Nh = N+altit ; 

return 2*tau*( N*Math.pow(c3*eccen/Nh, 2 . ) /(Nh* (1 .+eccen*tau) * (1 . -eccen*tau) ) -1. ) ; 

} 

/** Second derivative of T with respect to tau. 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn the second derivative of discrTO with respect to tau 
*/ 

double d2Tdtau2(f inal double tau, final double c3) 
{ 

final double N = curvN(tau) ; 
final double Nh = N+altit ; 
final double e2 = eccen*eccen ; 

final double Oneetau = (1 . +eccen*tau) * (1 . -eccen*tau) ; 

return -2.*( 1 . -N*c3*c3*e2* (1 . +3 . *altit*e2*tau*tau/Nh/0neetau) /Math. pow(Nh, 3 .) /Oneetau ) 

} 

/** Derivative d tau/d s [1/m] . 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn absolute value of the derivative according to eq (61) of the preprint. 
*/ 

double dtauds(double tau, double c3) 
{ 

return Math.sqrt( discrT(tau,c3) ) / (altit+curvM(tau) ) ; 

} 

/** Second derivative d~2 tau/d lambda~2. 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn value of the derivative of eq (64) of the preprint. 
*/ 

double d2taudlambda2 (double tau, double c3) 
{ 

/* The array contains values of E, sqrt(T) and l/(h+M) in the components facto[][0], 

* and their first derivatives with respect to tau in facto[][l]. 

*/ 

double [][] facto = new double [3] [2] ; 

facto[0][0] = GaussE(tau) ; 

facto[l][0] = Math.sqrt(discrT(tau,c3) ) ; 

facto [2] [0] = l./( curvM(tau)+altit) ; 

facto[0] [1] = dEdtau(tau) ; 

/* derivative of square root is one half of the inverse multiplied by interior derivative 
*/ 
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facto[l] [1] = dTdtau(tau,c3)/(2.*facto[l] [0] ) ; 

/* derivative of l/(h+M) is negative of its square multiplied by interior derivative 
*/ 

facto[2][l] = -dMdtau(tau)*facto[2] [0] *f acto [2] [0] ; 

/* The Leibniz rule of derivatives is applied to collect the derivative of 

* the product of f acto [i] [0] , i=0..2. 
*/ 

final double resul = f acto [0] [1] *f acto [1] [0] *f acto [2] [0] +f acto [0] [0] *f acto [1] [1] *f acto [2] [0] 
+f acto [0] [0] *f acto [1] [0] *f acto [2] [1] ; 

/* resul contains so far the derivative (d/dtau) (d tau/d lambda) , without 

* the c3 in the denominator. Postmultiply with d tau/ d lambda to generate 

* the d~2 tau/ d lambda"2. 
*/ 

return dtaudlambda(tau, c3) *resul/c3 ; 

} 

/** Third derivative d~3 tau/d lambda~3. 

* Spar am tau sine of the latitude 

* Oparam c3 the constant c3 along one individual trajectory 

* Oreturn value of the 2nd derivative of eq (64) of the preprint. 
*/ 

double d3taudlambda3 (double tau, double c3) 
{ 

final double T = discrT(tau,c3) ; 
/* Faa di Bruno formula: 

* d"3 tau/d lambda"3 = (d"2tau/d lambda~2)* (d/dtau) (dtau/dlambda)+ (dtau/dlambda) "2* (d~2/dtau~2) (dtau/dlambda 
*/ 

/* The array contains values of E, sqrt(T) and l/(h+M) in the components [] [0] , 

* and their first derivatives with respect to tau in [] [1] , 2nd derivatives in [] [2] . 
*/ 

double [][] facto = new double [3] [3] ; 

facto[0][0] = GaussE(tau) ; 
facto [1] [0] = Math.sqrt(T) ; 
facto[2][0] = l./( curvM(tau)+altit) ; 

facto[0][l] = dEdtau(tau) ; 

/* derivative of square root is one half of the inverse multiplied by interior derivative 
*/ 

final double Tprime = dTdtau(tau, c3) ; 
facto [1] [1] = Tprime/ (2. *f acto [1] [0] ) ; 

/* derivative of l/(h+M) is negative of its square multiplied by interior derivative 
*/ 

final double Mprime = dMdtau(tau) ; 
facto[2][l] = -Mprime*facto[2] [0] *f acto [2] [0] ; 

facto[0][2] = d2Edtau2(tau) ; 

/* First derivative of sqrt(T) was dTdtau/(2*sqrt (T) ) . The 2nd derivative therefore is 

* [ (d~2/dtau"2) T -(d T/d tau) "2/(2T)] / [2sqrt (T)] , which we store in facto [1] [2] . 
*/ 

facto[l][2] = ( d2Tdtau2(tau,c3) -Tprime*Tprime/( 2.*T ) ) /(2 . *f acto [1] [0] ) ; 

/* Second derivative of l/(h+M). First one was -(dM/dtau) /(h+M) "2 . Second one is 

/* - (d~2M/dtau~2) / (h+M) "2+2 (dM/dtau) "2/ (h+M) ~3 . 

*/ 

facto[2][2] = ( 2.*Mprime*Mprime*facto[2] [0] -d2Mdtau2(tau) ) *f acto [2] [0] *f acto [2] [0] ; 
/* The value of (d/d tau) (d tau/d lambda) . The implementation is not 

* tuned for quick evaluation, since the same value is calculated again in #d2taudlambda2 . 
*/ 

final double dtaudlambdadtau = ( f acto [0] [1] *f acto [1] [0] *f acto [2] [0] +f acto [0] [0] *f acto [1] [1] *f acto [2] [0] 
+facto[0] [0]*facto[l] [0] *f acto [2] [1] )/c3 ; 

/* the value of (d~2/dtau~2) (dtau/dlambda) with Leibniz rule, multinomial case 
*/ 
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final double d3taudlambdadtau2 = ( f acto [0] [2] *f acto [1] [0] *f acto [2] [0] +2 . *f acto [0] [1] *f acto [1] [1] *f acto [2] [0] 
+2 . *f acto [0] [1] *f acto [1] [0] *f acto [2] [1] +f acto [0] [0] *f acto [1] [2] *f acto [2] [0] 
+2 . *f acto [0] [0] *f acto [1] [1] *f acto [2] [1] +f acto [0] [0] *f acto [1] [0] *f acto [2] [2] ) /c3 ; 

/* Faa dl Bruno rule of derivatives is applied to collect the derivative 

* (d~2/dlambda~2) dtau/dlambda 

* = (d~2 tau/dlambda"2)*(d/dtau) (d tau/d lambda)+(d tau/d lambda) "2* (d"2/d tau"2) (d tau/d lambda) 
*/ 

return d2taudlambda2(tau, c3) *dtaudlambdadtau + Math.pow(dtaudlambda(tau,c3) ,2. )* d3taudlambdadtau2 ; 

} 

/** Second derivative d~2 s/d lambda"2 = (d/d lambda) (E/c3) . 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn value of the 2nd derivative of eq (63) of the preprint. 
*/ 

double d2sdlambda2 (double tau, double c3) 
{ 

/* First the derivative (d/dtau) (d s/d lambda) , without 

* the c3 in the denominator. Postmultiply with d tau/ d lambda to generate 

* the d~2 s/ d lambda"2. 
*/ 

return dtaudlambda(tau, c3) *dEdtau(tau) /c3 ; 

} 

/** Third derivative d~3 s /d lambda ~3. 

* This is currently implemented as a dummy to return 0. Effectively 

* the order of the approximation remains second order. 

* (Sparam tau sine of the latitude 

* (Sparam c3 the constant c3 along one individual trajectory 

* (Sreturn value of the third derivative of eq (67) of the preprint 
*/ 

double d3sdlambda3 (double tau, double c3) 
{ 

return 0. ; 

} 



/** Spherical approximation to parameter c3 

* (Sparam phi the latitudes at start and end [rad] 

* (Sparam lambd the longitudes at start and end [rad] 

* (Sreturn the parameter c3 obtained form the spherical approximation. 

* (Ssee eq (A19) in the preprint. 
*/ 

double c3Sphere (final doublet] phi, final double lambd[]) 
{ 

/* cosine of the angle Z */ 

final double cosZ = Math. sin(phi [0] ) *Math. sin (phi [1] ) 

+Math . cos (phi [0] ) *Math . cos (phi [1] ) *Math . cos (lambd [1] -lambd [0] ) ; 

return (altit+rho_equat) *Math. cos (phi [0] )*Math. cos (phi [1] )*Math. sin ( lambd [1] -lambd [0] ) /Math. sqrt (1 . -cosZ*cosZ) ; 

} 

/** Approximate sign d tau / ds 

* (Sparam phi the latitudes at start and end [rad] 

* (Sparam lambd the longitudes at start and end [rad] 

* (Sreturn the sign of the square root in eq (64) of the preprint, as obtained from 

* the spherical approximation. 
*/ 

double dtaudsSignum(f inal double [] phi, final double lambd []) 
{ 

/* cosine of the angle Z */ 

final double cosZ = Math. sin(phi [0] ) *Math. sin (phi [1] ) 

+Math . cos (phi [0] ) *Math .cos (phi [1] ) *Math . cos (lambd [1] -lambd [0] ) ; 
return Math. signum(Math. sin (phi [1] ) -Math. sin (phi [0] ) *cosZ) ; 

} 

/** Adjust longitude of the final point. 

* (Sparam lambd the longitudes at start and end [rad] 

* The value of lambd [1] is adjusted to the range lambd [0]+-pi to ensure that the 
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* trajectory chooses a path along the correct hemisphere. 
*/ 

static void ad j LambdaEnd (double lambd[]) 
{ 

final double dl = lambd [1] -lambd[0] ; 

if ( Math.abs( dl+2 . *Math.PI) < Math.abs(dl) ) 

lambd [1] += 2.*Math.PI ; 
if ( Math.abs( dl-2.*Math.PI) < Math.abs(dl) ) 

lambd [1] -= 2.*Math.PI ; 

} 

/** Nautical course in the local oblique horizontal [rad] 

* Oparam tau sine of the latitude 

* Oparam c3 the directional parameter of the solution to the inverse problem [m] 

* Oparam north sign of the square root, +1 or -1 

* Oreturn the direction angle (North over East) in the local oblique tangential plane [rad] 

* Osee eq (78) of the preprint 
*/ 

double naut Angle (final double tau, final double c3, final double north) 
{ 

/* sinkappa and coskappa are the sine and cosine of the angle multiplied by 

* a common positive factor, which retains the correct quadrant of the solution. 

*/ 

final double sinkappa = c3/(curvN(tau)+altit) ; 

final double coskappa = north* Math.sqrt( discrT(tau,c3) ) ; 

return Math. at an2( sinkappa, coskappa) ; 

} 



/** Compute one trajectory over the lambda interval, assuming c3 given. 

* Oparam phi start and (target) value of the geodetic latitude [rad] 

* Oparam lambd start and end value of the longitude [rad] 

* Oparam c3 fixed parameter of c3 

* Oparam Nsampl number of steps into which the interval lambd [0] -lambd [1] is divided 

* Oparam usampl if positive, each usampl 'th point of the result is printed 

* Oparam useRad one value of the AngleUnit to indicate which units are preferred for angles 

* Oparam taylOrd the order of the Taylor approximation in the FEM, must be 2 or 3. 

* Oreturn the overshooting of the value of tau versus the sine of phi[l] [rad] 
*/ 

double tauShoot (final double[] phi, final double lambd [] , final double c3, final int Nsampl, final int usampl, 
final int useRad, final int taylOrd) 

{ 

/* step width in the longitudes, including a sign */ 
final double dlam = (lambd [1] -lambd [0] ) /Nsampl ; 

/* current value of tau = sin(phi) , initialized with the value at the start */ 
double tau = Math. sin(phi [0] ) ; 

/* current value of s, the length along the line, initialized with the value at the start */ 
double s = 0. ; 

/* the sign which chooses a N or a S route. +1 or -1 */ 
double north = dtaudsSignum(phi , lambd) ; 

if ( usampl > ) 

System. out. println("\n\n") ; 

double kappa = naut Angle (tau, c3, north) ; 

System. out. println("# start course "+ kappa +" rad, " + Math. toDegrees (kappa) + " deg") ; 

/* loop over the finite elements */ 

for(int i=0; i < Nsampl; i++) 

{ 

/* in each of the elements, a 2nd order Taylor approximation 

* for dtau/ds is applied with (64) of the preprint 

*/ 

double derivl = north*dtaudlambda(tau, c3) ; 



/* below, the 'north' factor appears twice and cancels */ 
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double deriv2 = d2taudlambda2(tau,c3) ; 

double deriv3 = (taylOrd >= 3) ? north*d3taudlambda3(tau,c3) : 0.; 

/* in each of the elements, a 2nd order Taylor approximation 

* for ds/dlambda is applied with (64) of the preprint 
*/ 

double derivls = dsdlambda(tau, c3) ; 
double deriv2s = d2sdlambda2(tau, c3) ; 

double deriv3s = (taylOrd >= 3) ? d3sdlambda3(tau,c3) : 0.; 

/* Taylor extrapolation. Current value +(df /dx)*delta+(l/2) (d~2/dx~2)*delta~2+(l/6)*(d~3/dx~3)*delta~3 

* Update to the value at the right end of the interval. 
*/ 

tau += dlam*(derivl+dlam*(0.5*deriv2+dlam*deriv3/6.) ) ; 

/* We may have passed by a maximum of the equatorial distance 

* and sense this by comparing the location of the extremum on the current parabolic 

* estimator for tau(lambda) with the current interval. Setting the derivative of the 

* Taylor estimator to zero, (df /dx)+(d~2/dx~2)* (lambda-current lambda)=0, the 

* location of the extremum is -(df/dx)/( d~2f/dx~2) relative to the current lambda. 

* (Stodo implement the 3rd order Taylor Approximation as well. 
*/ 

final double lamFlat = -derivl/deriv2 ; 

/* Update of the path length along the trajectory in the Taylor approximation 
*/ 

s += dlam* (derivls+dlam* (0 . 5*deriv2s+dlam*deriv3s/6 . ) ) ; 

if ( usampl > ) 

if ( i 7, usampl ==0 I I i == Nsampl-1 ) 
{ 

/* longitude at right end of interval [rad] */ 
final double longi = lambd [0] +(i+l) *dlam ; 
double [] cart = getCartesian(tau, longi ) ; 

/* Print the x,y,z Cartesian coordinates as the first 3 columns per line */ 
System, out. print (" "+cart [0] + " "+cart[l] + " "+cart [2] ) ; 

kappa = naut Angle (tau, c3, north) ; 

/* Print longitude, geodetic latitude, nautical angle as columns 4 to 6 */ 
if ( useRad == AngleUnit .RAD ) 

System. out .print (" "+ longi +" "+ Math.asin(tau)+ " " +kappa ) ; 
else if ( useRad == AngleUnit .DEG ) 

System. out .print (" "+ Math. toDegrees (longi) + " " + Math. toDegrees (Math. asin(tau) ) 
+" "+Math. toDegrees (kappa) ) ; 

else 

System. out .print (" "+ 10 . *Math. toDegrees (longi) /9 . + " " 
+ 10 . *Math .toDegrees (Math . asin(tau) ) /9 . 
+" "+10. *Math. toDegrees (kappa) /9. ) ; 

/* Print length of trajectory in column 7 */ 
System. out. print (" "+s) ; 
System. out. printlnC") ; 

} 

if ( lamFlat/dlam >= 0. kk lamFlat/dlam < 1.) 
{ 

north *= -1. ; 

} 



/* tau now is the estimate of the final position; return the error */ 
return tau-Math. sin(phi [1] ) ; 

} 



/** Perform 4 iterations on adjusting the parameter c3. 

* Sparam phi the latitudes at start and end [rad] 

* Oparam lambd the longitudes at start and end [rad] 
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* (Sparam Nsampl the number of divisions of the longitude axis 

* (Sparam usampl undersampling factor for the printout in the converged trajectory 

* (Sparam useRad the unit of the angles which are reported 

* (Sparam taylOrd the order of the Taylor expansion of the core integral, 2 or 3. 

* (Sreturn a convergent to the parameter c3 
*/ 

public double c3shoot(f inal double[] phi, final double lambd[] .final int Nsampl, final int usampl, final int useRad, 
final int taylOrd) 

{ 

/* For a parabolic (2nd order) approximation of the shooting 

* error (units of tau) use 3 abscissa and target missing parameters 
*/ 

double [] err = new double [4] ; 
double [] c3spread = new double [4] ; 

/* start with the estimate from the spherical approximation, zero eccentricity */ 
c3spread[0] = c3Sphere(phi , lambd) ; 

err[0] = tauShoot (phi , lambd, c3spread [0] , Nsampl , -1 , useRad, taylOrd) ; 

System. out .println("# c3 (spherical) "+c3spread [0] +" err "+err[0]) ; 

/* A stepping parameter for the variation */ 
final double c3Delt = -c3spread [0] *0 . 005 ; 

c3spread[l] = c3spread [0] +c3Delt ; 

err[l] = tauShoot (phi , lambd, c3spread [1] , Nsampl , -1 , useRad, taylOrd) ; 
/* estimate of the 3rd c3 by a linear interpolation 

* between the previous 2 to get zero error. The ansatz is err=err [0] + (err [1] -err [0] ) /c3Delt* (c-c [0] ) = 0, 

* where (err [1] -err [0] ) /c3Delt is the slope of the error curve, which is solved for c. 
*/ 

c3spread[2] = c3spread [0] -c3Delt*err [0] / (err [1] -err [0] ) ; 

err[2] = tauShoot (phi , lambd, c3spread [2] , Nsampl , -1 , useRad, taylOrd) ; 

System. out. println("# c3 (linear extrapol) "+c3spread [2] +" err "+err[2]) ; 

/* Lagrange quadratic interpolation between (c3spred[0] ,err [0] ) , (c3spread [1] , err [1] ) 

* and (c3srepad[2] ,err [2] ) : err = eof Cinter [0] +eof Cinter [1] *c3spread+eof Cinter [2] *c3spread~2 

* up to a common factor which is omitted. 
*/ 

double [] eof Cinter = new double [3] ; 

eof Cinter [0] = err [2] *c3spread [0] *c3spread [1] * (c3spread [0] -c3spread [1] ) 

+err [1] *c3spread[0] *c3spread[2] * (c3spread[2] -c3spread[0] ) 

+err [0] *c3spread[l] *c3spread[2] * (c3spread[l] -c3spread[2] ) ; 
eof Cinter [1] = err [2] * (Math.pow(c3spread [1] , 2 . ) -Math.pow(c3spread [0] ,2.)) 

+ err[l]*(Math.pow(c3spread[0] , 2 . ) -Math.pow(c3spread [2] ,2.)) 

+ err[0]*(Math.pow(c3spread[2] , 2 . ) -Math.pow(c3spread [1] ,2.)) ; 
eof Cinter [2] = err [2] * (c3spread [0] -c3spread [1] ) 

+err [1] * (c3spread [2] -c3spread [0] ) 

+err [0] * (c3spread [1] -c3spread [2] ) ; 



/* Set this quadratic to zero; normalize to monic polynomial */ 
eof Cinter [0] /= eof Cinter [2] ; 
eof Cinter [1] /= eof Cinter [2] ; 

eofCinter[2] = Math. sqrt (0 . 25*Math.pow(eof Cinter [1] , 2) -eof Cinter [0] ) ; 
/* From the two solutions take the one closer to the previous estimate */ 

c3spread[3] = (c3spread [2] +0 . 5*eof Cinter [1] > 0.) ? -0 . 5*eof Cinter [1] +eof Cinter [2] : -0 . 5*eof Cinter [1] -eof Cinter [2] ; 
err[3] = tauShoot (phi , lambd, c3spread [3] , Nsampl , usampl , useRad, taylOrd) ; 
System. out .println("# c3 (quadratic extrapol) "+c3spread [3] +" err "+err[3]) ; 

return c3spread[3] ; 

} 

/** 

* usage : 

* <tt> 

* java Geod [-h altitude/m] [-s #longitude steps] [-u sfac] [-R equat radius/m] [-e eccent] [-r] 

* [-g] [-2] phil laml phi2 lam2 
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* </tt> 

* <br> 

* Command line parameters 

* <ul> 

* <li> 

* -h followed by a floating point number specifies the altitude of the 

* trajectory above the ellipsoid. This is given in the same units as those of the 

* parameter of the equatorial radius. The default is zero. 

* </li> 

* <li> 

* -s followed by an integer denotes the number of steps between the starting and 

* final point of the trajectory in the finite element method. The default is 400. 

* </li> 

* <li> 

* -u followed by an integer specifies an undersampling factor for listing the 

* parameters of the trajectory on standard output. The elements actually listed 

* are given by the sampling steps divided by the undersampling factor. 

* The default is 10. 

* </li> 

* <li> 

* -R followed by a floating point number defines the equatorial radius of 

* the surface of the ellipsoid in the same units as the option for the altitude. 

* The default is the WGS84 value. 

* </li> 

* <li> 

* -e followed by a floating point number defines the eccentricity of the 

* surface of the ellipsoid. The default is the WGS84 value. 

* </li> 

* <li> 

* -r without any additional number specifies that the four mandatory command 

* line parameters are provided in units of radians. The default is degrees. 

* </li> 

* <li> 

* -g without any additional number specifies that the four mandatory command 

* line parameters are provided in units of gons. The default is degrees. 

* </li> 

* </ul> 

* There are four mandatory parameters, the goedetic latitude and longitude 

* of the start of the trajectory, and the geodetic latitude and longitude of 

* the final point of the trajectory. 

* <br> 

* The standard output contains comment lines (starting with hashes) and 

* one line per point on the trajectory with the following columns (separated by blanks): 

* <ul> 

* <li> 

* The first three colums are the x, y and z Cartesian coordinate of the point in the 

* same units as with the equatorial radius and altitude parameters. 

* </li> 

* <li> 

* Columns 4 and 5 are the longitude and geodetic latitude of the point 

* in the same units (deg or rad) as the command line parameters. 

* </li> 

* <li> 

* Column 6 is the nautical angle in the local topocentric tangential plane 

* in the same units (deg or rad or gon) as the command line parameters. 

* </li> 

* <li> 

* Column 7 is the distance along the geodetic 

* in the same units as the command line parameters. 

* </li> 

* </ul> 

* <br> 

* Examples 

* <ul> 

* <li> 

* <tt> 

* java Geod -s 1000 -e 0.08227185417541244347812117 -R 6378206.4 8.973611111 -79.573333333333 21.435000000 -158.02583333 

* </tt> 

* computes the trajectory for the case from Panama to Hawaii by Paul D. Thomas, 

* J. Geophys. Res. vol 70 (1965) p 3331 in his Figure 3. 
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* <li> 

* <tt> 

* java Geod -R 637838.799919243 -g -e 0.0815815883368028 54.262654 11.977469 115.959876 

* </tt> 

* computes the trajectory for the case from Paris to Saigon by H. M. Dufour 

* in Bull. Geo. vol 32 (1953) p 26, with f =1/300, angles in gons . 

* </ul> 

* Osince 2008-04-12 
*/ 

public static void main(String[] args) throws Exception 
{ 

/* default equatorial radius set to WGS84 */ 
double rho = WGS84_RH0_E ; 

/* default eccentricity set to WGS84 

* that is 0.08195646813308029492885488 ; 

*/ 

double e = (2.-1. /WGS84_FLAT) /WGS84_FLAT ; 

/* default altitude above ellipsoid is zero */ 
double h = 0. ; 

/* default number of elements */ 
int Wsampl = 400 ; 

/* Default undersampling factor for the printed output. 
*/ 

int usampl = 10 ; 

/* Default Taylor order in the FEM is 3. 
*/ 

int taylOrd = 3 ; 

/* no defaults for the coordinates */ 
double [] phi = new double [2] ; 
double [] lambd = new double [2] ; 

/* units of final angular arguments. */ 
int useRad = AngleUnit .DEG ; 

/* parse options and parameters 
*/ 

for(int argc =0 ; argc < args. length; argc++) 
{ 

if ( args [argc] . compareTo("-h") == 0) 

h = Double. parseDouble (args [++argc] ) ; 
else if ( args [argc] . compareTo("-s") == ) 

Nsampl = Integer .parselnt (args [++argc] ) ; 
else if ( args [argc] . compareTo("-u" ) == ) 

usampl = Integer .parselnt (args [++argc] ) ; 
else if ( args [argc] . compareTo("-R" ) ==0 ) 

rho = Double. parseDouble (args [++argc] ) ; 
else if ( args[argc] .compareTo("-e") ==0 ) 

e = Double. parseDouble (args [++argc] ) ; 
else if ( args[argc] . compareTo("-r" ) ==0 ) 

useRad = AngleUnit .RAD ; 
else if ( args [argc] . compareTo("-g" ) ==0 ) 

useRad = AngleUnit . G0N ; 
else if ( args [argc] . compareTo("-2" ) == ) 

taylOrd = 2 ; 

else 
{ 

phi[0] = Double. parseDouble (args [argc++] ) ; 
lambd [0] = Double .parseDouble (args [argc++] ) ; 
phi[l] = Double. parseDouble (args [argc++] ) ; 
lambd [1] = Double .parseDouble (args [argc++] ) ; 

} 

} 



/* If user provided degrees or gons on the command line, convert to radians for further processing 
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* 400 gons are 360 degrees. 
*/ 

switch( useRad) 
{ 

case AngleUnit .DEG : 

for(int i=0 ; i < 2 ; i++) 
{ 

phi [i] = Math. toRadians (phi [i] ) ; 
lambd[i] = Math. toRadians (lambd [i] ) ; 

} 

break ; 
case AngleUnit .GQN : 

for(int i=0 ; i < 2 ; i++) 
{ 

phi[i] = Math. toRadians (360. *phi[i]/400.) ; 
lambd [i] = Math. toRadians (360 . *lambd[i] /400 . ) ; 

} 

break ; 

} 

/* normalize input such that the final longitude is on the same hemisphere */ 
adjLambdaEnd(lambd) ; 

/* log the input parameters */ 

System. out .println("# equat Radius "+ rho+" , eccentricity "+e+", altitude "+h) ; 
System. out .println( "# inverse flattening "+ 1 . /(I . -Math . sqrt (1 . -e*e) ) ) ; 

System. out .println("# start "+phi[0] + " rad, "+lambd[0]+" rad; end "+phi[l] + " rad, "+lambd[l]+" rad") ; 

switch( useRad) 

{ 

case AngleUnit .DEG : 

System. out. println("# start "+Math. toDegrees (phi [0] ) + " deg, "+ Math. toDegrees (lambd [0] ) 

+" deg; end "+Math.toDegrees(phi [1] ) + " deg, "+Math. toDegrees (lambd [1] )+" deg") ; 
break ; 
case AngleUnit .GON : 

System. out. println("# start "+10 . *Math. toDegrees (phi [0] )/9 . + " gon, "+ 10 . *Math. toDegrees (lambd [0] )/9 . 
+ " gon; end " + 10 . *Math. toDegrees(phi [1] ) /9 . + " gon, "+10 . *Math. toDegrees (lambd [1] )/9 . + " gon") ; 

} 

System. out .println("# "+Nsampl + " elements, Taylor order " + taylQrd) ; 

/* define the shape of the ellipsoid */ 
Geod g = new Geod(rho,e,h) ; 

/* Determine the main parameter of the bundle of geodesies through initial point to 

* solve the inverse problem, and print parameters along the trajectory 

* in the 4th iteration. 
*/ 

g. c3shoot (phi , lambd, Nsampl ,usampl , useRad, taylOrd) ; 

} 

} /* Geoid */ 
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